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Abstract 

In molecular reactions at the microscopic level the appearance of resonances 
has an important influence on the reactivity. It is important to predict when a 
bound state transitions into a resonance and how these transitions depend on 
various system parameters such as internuclear distances. The dynamics of such 
systems are described by the time-independent Schrodinger equation and the 
resonances are modeled by poles of the .S'-matrix. 

Using numerical continuation methods and bifurcation theory, techniques 
which find their roots in the study of dynamical systems, we are able to develop 
efficient and robust methods to study the transitions of bound states into 
resonances. By applying Keller's Pseudo-Arclength continuation, we can minimize 
the numerical complexity of our algorithm. As continuation methods generally 
assume smooth and well-behaving functions and the S'-matrix is neither, special 
care has been taken to ensure accurate results. 

We have successfully applied our approach in a number of model problems 
involving the radial Schrodinger equation. 
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1. Introduction 

Over the last couple of decades several reliable numerical methods have been 
developed for continuation of solutions and bifurcation analysis for dynamical 
systems [1, 2, 3, 4]. In this contribution we investigate the application of these 
methods in the context of quantum physics. In particular, we use numerical 
continuation to trace the dependence of the energy and width of resonances on 
the system parameters. This is relevant e.g. in low energy electron-molecule 
scattering where the occurrence and structure of the resonance depend on the 
internuclear distance in the molecule. 

Our system of interest fits the radial Schrodinger equation, a subclass of 
Sturm-Liouville boundary value problems. For these types of problems, there 
exist several very accurate methods that find the bound state eigenvalues [5, 6]. 
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In many physical systems, however, it is also valuable for finding the resonant 
states that have a complex valued energy. 

We define resonances and bound states as solutions of the Schrodinger 
equation for an energy where the S-matrix has a pole [7, 8, 9]. The S-matrix is a 
function of the complex momentum k that can be extracted from the solution at 
the end of the domain. It also depends on the system parameters. We introduce 
a regularization procedure that transforms the poles into zeros and smoothes 
the behavior near the origin in the fc-plane. This allows the application of the 
pseudo-arclength continuation method to trace the trajectory of the zeros, and 
hence the poles, as the system parameter changes. 

The outline of the paper is as follows. In section 2 we present an overview 
of the concepts underlying the numerical method that constructs a solution set 
of a non-linear equation with the help of numerical continuation. As indicated 
in the application in section 4, we use an implementation of these methods 
provided by the AUTO package [10]. In section 3 we review the concepts related 
to the Schrodinger equation, its solution through the renormalized Numerov 
method and the extraction of the S'-matrix from the numerical wave function. It 
is the poles of the S'-matrix that are subjected to the numerical continuation 
methods of section 2. Finally, in section 4 we demonstrate our approach on two 
models describing a single-particle in three dimensions in a spherically-symmetric 
potential. Using partial wave expansion, these scattering problems reduce to a 
radial Schrodinger equation. 

2. Numerical continuation methods 

Numerical continuation methods approximate the solution set of some non- 
linear equation -F(u, A) = that depends on a system parameter A: 

F:R n+1 — >R n :(u,A) i — >F(u,\), (1) 

where u £ M. n . The Implicit Function Theorem states that under certain 
continuity conditions the solution set is a one-dimensional manifold and can be 
parameterized by some real parameter s. The choice of that parameter is an 
important one and depends on the method used. Generally, we are interested in 
the evolution of the solutions u in terms of A and this suggests to take A as the 
continuation parameter. However, this choice may result in difficulties when the 
solution path passes through a fold. The pseudo-arclength continuation [1] deals 
with these situations gracefully. 

We introduce several notations used throughout this paper. When the 
distinction between the function variables u and the parameter A is irrelevant, we 
write x = (u, A) and Xj = (uj, Ai) for the subsequent points on the solution curve. 
The continuation curve is denoted by x(s), which emphasizes the dependence 
on the continuation parameter s. The initial point on the curve is associated 
with s — and written as Xo = x(0). Numerical continuation methods use this 
point on the curve, along with an initial direction of continuation to construct a 
sequence of points 

{ Xi | i = 0,...,N and F(x,-) = 0} , (2) 
that approximates the solution curve. 
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Figure 1: Graphical representation of one predictor-corrector step in pseudo-arclength continu- 
ation as discussed in section 2. 



2.1. Pseudo-Arclength Continuation 

The algorithm follows a predictor- corrector scheme to construct, starting 
from an initial solution point Xo, the successive points on the solution curve. 

The predictor step is an Euler predictor that uses the unit length tangent 
vector Xj to the curve at a solution point x^ (thus satisfying F(xi) = 0) and a 
step size As to predict a guess x^ +1 for the next point on the curve: 

xf+i = x, + Asx. t . (3) 

The corrector step improves the guess x^ +1 with a Newton iteration on the 
augmented system to obtain a new solution point Xj+i. This augmented system 
has, in addition to the constraint .F(xj+i) = 0, the requirement that Xj+i must 
lie on the hyperplane through x^ +1 perpendicular to Xj, the tangent to the 
previous solution. This translates to 

= (4) 
(x l+ i - xf +1 ) • Xi = 0. 

This system is a map from to E" +1 and defines, under some conditions that 
are usually met, uniquely the next point on the solution curve. It is the point 
of intersection between the hyperplane and the curve shown in figure 1. These 
steps are common to other Euler-Newton like methods and other approaches to 
define the next point on the curve are discussed in [11]. 

The tangent vector for the next step is computed by solving: 



F u (x i+ i) F\(x i+ i) \ / Ui+i \ _ f 



(5) 



and normalizing ||xj_|_i|| = 1. The right-hand side of equation (5) is a column 
vector consisting of zeros except on the last row. 

Note that the Jacobian, F x (xj+i), is required both for the calculation of the 
tangent direction and for the calculation of the Newton corrections. In our 
application we only have F numerically so we need to approximate the Jacobian. 
This is done using finite differences. The jth column of the Jacobian matrix is 
found by a central difference and requires two solutions with slightly different 
arguments: 

, ( v» F(x i+1 + ee 3 -) - F(x i+ i - eey) 

[*x(Xi+i))j ^ ) \P) 
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where is is the jth unit vector. A discussion on the optimal choice of e given 
the machine precision is found in [12]. 

2.2. Regular and singular solutions 

An important notion is the regularity of a solution point [4, 3]: a point 
Xj £ M. n+1 on the solution curve is a regular solution of F(x) = if the Jacobian 
matrix, f x (xj), has maximal rank. Otherwise the solution is singular. Since 
F x (x.i) has n rows and n + 1 columns, its maximal rank is n. 

Another important notion is bifurcation. The solution is said to bifurcate [13] 
from the solution u t at a parameter value A t if there are two or more distinct 
solutions which approach u t as A tends to a threshold value A(. A more rigorous 
definition of a bifurcation point can be found in [14]. 

The connection between these two definitions is that a bifurcation point x, 
of -F(x) = must be a singular solution which means that: 

rank(F x (x i )) < n, (7) 

and consequently (rank- nullity theorem): 

dimkcr(F x (x l )) > 2. (8) 

In case the equality in (8) holds, we call Xj a simple bifurcation point [14]. We 
assume this is the only type of bifurcation that occurs in the systems we study 
here. 

Following [14] we detect these bifurcation points by looking at the sign of 
the determinant of the augmented Jacobian matrix. When traversing a solution 
branch a simple bifurcation point lies between two solutions and Xj+i if and 
only if 

sgndet ( ) ^ sgndet ( F ^+ l] ) . (9) 

This allows to find the bifurcation point accurately with a straightforward yet 
rather slow convergence procedure. Note that Xj and Xj+i must be close to each 
other to avoid "overshooting" bifurcation points. 

2. 3. Branching 

When two solution curves meet in a simple bifurcation point x t , the dimension 
of the nullspace of F x (x t ) is two. This nullspace is then spanned by two 
orthonormal vectors ti and t 2 . At the same time, the left nullspace of F x (x t ) 
is one dimensional since F is a function from W l+1 to E™. It is spanned by a 
vector n\. 

The two tangent vectors to the curves that depart from the bifurcation point 
can now be written as a linear combination x t = at\ + fiti of the vectors that 
span the nullspace. Since the curves fit F(x(s)) = 0, we can differentiate twice 
to s and find that the tangent directions fit 

F xx x t x t + F x x t = 0. (10) 

Projection on n\ leads to the algebraic bifurcation equation [3, 15]: 

Cna 2 + 2C 12 af3 + C 22 f3 2 = 0, (11) 
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with Cn = nj F xx titi, C\2 = n\ F xx tit2 and C22 — F xx t2t2- In addition we 
have a 2 + 1 = 1, since the tangent vectors are normalized. 

The construction of this equation requires a numerical calculation of the 
Hessian in the bifurcation point and the determination of the vectors that span 
the nullspaces. The solution of the algebraic bifurcation equation gives us the 
tangent vectors to the curves that depart from the bifurcation point. 

2.4- Implementation 

During initial prototyping we have implemented the above methods in Mat- 
LAB. For the development of production code we have relied on the well-known 
implementation of these algorithms in the AUTO package [10, 16]. An alternative 
implementation is available in the LOCA package which is part of the Trilinos 
project [17]. 

3. Quantum scattering concepts 

In this section we review some of the concepts related to the solution of the 
time-independent Schrodinger equation through partial wave analysis and to the 
S'-matrix and its properties. 

3.1. The radial Schrodinger equation and the S -matrix 
The time-independent Schrodinger equation 

f-~A + V(f,A)) iP(r) = EiP(r), (12) 

describes the states of a quantum system with potential V at energy E. We 
let the potential depend on a parameter A. How the potential depends on the 
parameter A is arbitrary. Any choice is acceptable provided the A-dependence 
is smooth. One possible choice is to scale the potential with a strength A as in 
V(f,A) = XV (f). 

In almost all physically relevant situations, V is spherically symmetric, i.e. 
a function of the radial coordinate r = \r\ only. One then transforms equation 
(12) to spherical coordinates (r, 9, ip) and applies the method of separation of 
variables — partial wave analysis in physics parlance — to solve as [18, 19]: 

i>(r) = J2 c irn^Y lm (9,<p), (13) 

The spherical harmonics Yi m are the solutions to the angular equation that is 
independent of V. The integer I is the angular momentum. For each I the is 
determined by a radial equation of the following form: 

H I* + y(r ' A) + ^ (r) = EMr) > (14) 

One refers to the sum of V and the /-dependent term as the effective potential. 
This equation belongs to a subclass of Sturm-Liouville boundary value problems 
with p(x) — 1, w(x) — 1 and q(x) equal to the effective potential. 
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We assume that for V(r, A) the following holds 

V(r, A) = 0(r~ 3 ~ e ) for r ->• oo and e > , , 

V(r, A) = 0(r- 2+c ) for r ^ and e > 0, ^ ^ 

V decays faster than r -3 at infinity and is less singular than r~ 2 at the origin. The 
requirement at infinity limits us to so-called short-range potentials. Extending our 
approach to the class of long-range interactions requires substantial modifications 
and is an important direction for future work. 

The solution -0/ of (14) needs to fit the homogeneous Dirichlet boundary 
condition at r = 0. Because of the short range of the potential, the solution 
becomes at r oo a linear combination of the two fundamental solutions of the 
free Schrodingcr equation (i.e. without potential term V) 

+ -*)*(-)-* < 16 > 

where k — \f2E £ C is the complex momentum. The fundamental solutions are 
the spherical Riccati-Bessel and Riccati- Neumann functions [20, 7]. Thus, in the 
asymptotic region we have for the solution of (14) 

ipi(r)^-Ai(k,X)hf(hr) + B l (k,X)hl(kr) for r -> oo, (17) 

with Ai(k, A), Bi(k, A) € C. These constants depend on the momentum k and 
system parameter A. The hf are Riccati-Hankel functions of the first (+) and 
second (— ) kinds. These functions behave asymptotically (up to a phase) as 
e lkr , an outgoing, and e~ lkr 1 an incoming wave. The solution is then interpreted 
as a superposition of an incoming (hf) and an outgoing (hf) wave. 

As the solution is only defined up to an overall normalization, we can 
renormalize it as follows [7]: 

Mr) = l(hT(kr)-S l (k,X)h+(kr)), with Si(k,X) = (18) 

This introduces the ^-matrix, a function of the momentum k and depending on 
A. It determines the phase of the outgoing, scattered wave w.r.t. the incoming 
wave. 

3.2. Resonances and bound states as poles of the S -matrix 

It is well established and discussed in several textbooks that poles of Si(k, A) 
with 3?(fc) > correspond to bound states with energy fc 2 /2 and poles with 
^s(k) < correspond to resonances. 

Indeed, if E is a negative real number where Si diverges, it means that the 
solution is asymptotically a multiple of hf only. Since E < 0, the momentum k 
is purely imaginary and ht becomes a decaying exponential that fits the zero 
boundary conditions as r — > oo. The solution then fits the boundary conditions 
V>(0) = and ip(oo) = and is then a bound state solution of the boundary 
value problem. 

On the other hand, if E is a complex number with Q(k) < where Si has a 
pole, the state is classified as a resonance. Again the asymptotic solution is a 
multiple of hf, an outgoing oscillating wave, only. 
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Figure 2: Schematic trajectories of the poles of the S-matrix for the Gaussian well example 
with 1 = 1. The four parts of the solution meet in the origin of the complex fc-plane for 
A ~ 6.0497. In this article this point is interpreted as a bifurcation point. 



Furthermore, as the system parameters change, the poles can move from the 
upper part towards the lower of the complex plane along a continuous curve 
and the solutions then transition from bound to resonant state. A thorough 
discussion on such trajectories of poles is given in [7, 8, 9]. 

In mathematical terms, we are faced with the problem of finding fc(A) : K — > C, 
the solution set of 

5 i (fc,A)- 1 =0. (19) 

When we treat its real and imaginary parts as independent variables, Si(k, A) is a 
function from M 3 to K 2 . It is a non-linear function and only for a few potentials 
it is available as an analytical expression. In general, the value of the S-matrix 
for a given 5ft(fc), Q(k) and A is only found through the numerical solution of the 
Schrodinger equation with methods such as the i?-matrix [21], J-matrix [22], 
ECS [23], shooting methods [24] and others. Given a numerical expression for 
the wave function, we extract the Ai(k, A) and Bi(k, A) from the solution with 
the help of the Wronskian W at the asymptotic boundary R of the domain. 
More specifically, from (17) we derive: 

wUm, h-(kR)) 

Mk, A) - — A- (20) 

W[h+{kR), h^(kR)j 



Bi(k, A) = A- ; \. (21) 

(h+(kR), h-(kR) 



w 

This then leads to the following expression for the S'-matrix 



Si(k, A) = — ) f , (22) 

W[MR), h+(kR)j 

that can be computed numerically, provided the first derivative of the wave 
function at R, the asymptotic boundary of the domain, can be computed. 

Note that the different solution curves of (19) can meet each other in a 
single point. For a radial equation with I = 1, for example, the situation is well 



understood and we illustrate this with the help of figure 2. As A — > At, the 
critical system parameter where a bound state becomes a resonance, the pole 
moves down on the positive imaginary axis towards the origin. At the same time, 
another pole corresponding to the virtual state, approaches the origin along 
the negative imaginary axis from below. At A t the two poles of the S-matrix 
coalesce into a single double pole in the origin. For system parameters beyond 
At there are again two separate poles corresponding to resonances. They lie in 
the third and fourth quadrant of the complex fc-plane. 

3. 3. Transforming poles into zeros 

The numerical method discussed in section 2 requires some smoothness 
conditions on the function -F(x) to work in a reliable and fast way. To achieve a 
quadratic convergence rate during the Newton correction it is well known that 
the Jacobian F x needs to be locally Lipschitz continuous. 

We intend to apply the continuation method to track the path of the zeros 
of Si(k, A)" 1 as A varies. Unfortunately, this function is meromorphic for the 
potentials of interest [7] and does not fit these smoothness requirements, especially 
for \k\ <C 1. Indeed, the scattering matrix has the property that Sj(fc*,A) = 
Si(k, A) -1 , where k* is the complex conjugate of k. This means that if Si(k, A) -1 
has a zero in some fco, it will also have a pole in k$. And as A approaches At, fco 
moves towards the origin. In this situation a zero in fco and a pole in fcp approach 
each other and at the critical parameter At they will coalesce. It is clear that 
in a neighbourhood around the critical point, |fco| <C 1 and |A — At] *C 1, the 
derivatives of S^fcjA) -1 cannot satisfy these smoothness conditions. 

In order to desingularize Si(k, A) -1 , i.e. to avoid this deteriorating behavior 
as |fc| — > 0, we transform to a new function, related to the 5-matrix, but having 
polynomial behavior for |fc| — » 0. This function is defined as 

1.21 + 1 

with I the angular momentum. It is clear that for fc ^ 0, Fi will have a zero if 
and only if Si has a pole. 

Furthermore, we can show that this function is proportional to the Jost 
function J-)(fc, A), familiar from scattering theory [7, 25, 26]. It is related to the 
ratio of the regular solution tpi(k, r) and the normalized solution ipi(k,r) and 
which is an analytic function for a wide class of potentials and behaves as a 
polynomial around the origin of the complex plane. To show this proportionality 
we use equations (11.19) and (12.145) from [8]. We have that 



Si(k 1 \) = l-Aik- 1 drkrji(kr)V(r,\)i>i(k,r), (24) 
Jo 

where ipi(k, r) is the solution to equation (14). This wave function is proportional 
to the regular solution ipi(k, r) 

Mk > r)= Fdk,\)(2l + !)!! ■ (25) 
where !! indicates the double factorial [19]. 



8 



For I = this ipi(k,r) fits equation (14) with boundary f(k,0) = and 
ip'(k, 0) = 1. In the case I > the condition is: lim r _j.o r~ l ~ 1 (pi(k, r) = 1. 
With the help of [8], it is clear that 

fi(M) = ^(MM = ~ (26) 

where 

C = rrrTTTw , ^„ / dr fcr j,(fcr)7(r, A)w(*,r). (27) 



fc'+ 1 (2/ + l)!! 7 

This constant is bounded. Indeed, we have the bound from [7] and [27] 

\krMkr)\<C^^y\^K (28) 
and in a similar way from [8] we have 

lw(*,r)| < C 2 ( e MV\r\k\-i-\ (29) 



(30) 



1 + \k\r 

with constants C\ and C%. So we get 

21+2 



This bound is finite if the integral over the potential is finite. As we can see 
from (27) it is clear that as k — > 0, C is only zero for very specific potentials. 
We conclude that Fi(k,r) is bounded for a wide range of problems. 

Note that tracking the zeros of Ai(k, A) or Bi(k, A) is not an alternative since 
these functions also suffer from the presence of poles. These poles are removed 
by taking the S- matrix, the ratio of the Ai and B%. 

3.4. The tangent directions in the bifurcation point 

An advantage of working in the fc-plane, instead of the -E-plane, is that the 
tangents to the solutions that emerge from the bifurcation point are orthogonal 
for problems with I > 1. 

Indeed, around k = 0, the Jost function can be expanded in the form [8] 

Fl{k, A) = ai(A) + a 2 (A)fc 2 + . . . + &!(A)£; 2i+1 + 6 2 (A)fc 2 '+ 3 + . . . (31) 

where the coefficients are real functions of the system parameter. Around 
(0, 0, At), with I > 1 this can further written as 

F{k, A) = a(X - A t ) + /3fc 2 + 0(k s ) +0((X- A t ) 2 ) , (32) 

where a, f3 € K.. When we write k = x + iy and define the function 

(x,y,X) ^ (n(r(x + iy,\)) t Qt(r(x + iy,\))), (33) 
this equation G(x, y, A) = follows the full problem up to order fc 3 and (A — A t ) 2 . 
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The Jacobian in the point (x, y, A) is then 



which obviously reduces to the following rank one matrix at the bifurcation point 
x t = (0,0, A t ) 

«•-(!! US)- (35) 

A basis for the ker (G x (x t )) is then 

ker(G x (x t )) = span | ( j , f 1 j | , (36) 
and for ker(G x ) it is 

ker (G x (x f ) T ) = span {(?)}■ ( 37 ) 

The Hessians are then 

Gxx = ( ( 2 ) ' (2 2 ) ' (0 ) ) (38) 

The coefhcients of the algebraic bifurcation equation (11) are then C\\ = 0, 
C12 = 2, C22 = what leads to the equation to be solved: 

a/3 = and a 2 + /3 2 = 1. (39) 

The solutions, which are (a,f3) = (0,1) or (1,0) then lead to two tangent 
directions at the bifurcation point: x^i = (1,0, 0) T and x t 2 = (0,1, 0) T . The 
direction x t i corresponds to the two resonances that leave along the real axis 
and x t 2 is the direction along the imaginary axis from which the bound and 
anti-bound states approach the threshold. These tangent vectors are indeed 
orthogonal. 

If we would use numerical continuation in the i?-plane, these two tangent 
directions would coincide. 

Note that AUTO solves, when it detects a bifurcation point, the algebraic 
bifurcation equation numerically. 



4. Numerical application 

4-1. Implementation 

For testing purposes we have developed an implementation of the algorithm 
described. It consists of two main parts: 

1. A solver for the Schrodinger equation and the associated routines to obtain 
a numerical approximation of the S'-matrix. As indicated in section 3.2 
many solvers are possible, each suitable for a range of potentials, domains 
or dimensionality of the problem. For the two examples we present, dealing 
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with the radial Schrodinger equation, the Numerov method has proven 
very successful. 

The Numerov method [28, 24] is a shooting method that exploits the 
absence of first order terms in the Schrodinger equation to arrive at a fairly 
straightforward algorithm, with equidistant steps, that is of fourth order. 
In [24] the renormalized Numerov method was proposed, a reformulation 
of the algorithm in terms of the ratio of the wave function in successive 
grid points. We have implemented this algorithm in C++ code. 
The first derivative of the wave function, required to compute the ratio 
of the Wronskians in (22) is determined with a formula given in [24] and 
retains the order 0(h A ). We have tested this convergence behavior in 
our implementation and found that it holds except for potentials with 
discontinuities, such as the square well. However, this does not prevent 
the application of the method: it simply lowers the convergence rate. 
2. A routine that performs the numerical continuation process with detection 
of branches. For this purpose, we use the well-known library for numer- 
ical continuation AUTO [16]. The numerical routines for the necessary 
computation of the Jacobian matrix are also provided by AUTO and use a 
second order central difference scheme. A comparative study with other 
continuation libraries is under consideration. 

4-2. Gaussian potential 

As a first model problem we take the third partial wave (I = 3) in a Gaussian 
potential well: 

V(r,X) = ~Xe- r \ (40) 

For potential strength A = 188 the system has 6 bound states. Decreasing the 
potential strength pushes these bound state energies towards zero and transforms 
them successively into resonances. For the renormalized Numerov solver an 
integration grid r e [0, 4.8] with 8192 points was used. At the end of this interval, 
the influence of the Gaussian potential is smaller than 10~ 7 and is considered 
negligible. The shooting method is started with the boundary condition at r = 0, 

iMo) = o. 

The starting points x = (fcsjf,fccj,A) for the six branches were chosen on 
the positive imaginary axis of the fc-plane, in a region close to the origin to 
ensure convergence of the solver. They are presented in table 1. To confirm our 
renormalized Numerov values we have also computed them using the CPM{ 16,14} 
method implemented in matslise [6]. There are no significant differences. 

Continuation was started in these points in the direction of the origin with 
an initial prediction step As = 1CP 2 which may vary dynamically between 10 -4 
and 5 x 10 -2 . The critical transition points, where the continuation branches 
off, were found at the origin of the fc-plane for threshold values for A given in 
table 2. 

The resulting trajectories of the continuation process are shown in figures 
3(a), 3(b) and 4. The time to compute each trajectory is of the order of several 
seconds on modern desktop computer hardware, depending on the step size and 
the number of continuation points. 
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n 


A 




Zs(k)(Ren.Numerov) 




S(fc) (matslise) 





25 


9 


343034507158935e-01 


9 


343034516458660e-01 


1 


46 


1 


226422927658922e+00 


1 


226422927676387e+00 


2 


72 


1 


207656897946988e+00 


1 


207656897794478e+00 


3 


104 


1 


174028026341686e+00 


1 


174028025751143e+00 


4 


142 


1 


125495438561443e+00 


1 


125495437195381e+00 


5 


188 


1 


294921256799873e+00 


1 


294921252331416e+00 



Table 1: Starting points for the six continuation branches in the Gaussian well example. 
3t(fc) = for all points. 



n 




A 




SR(fc) 




3 


(fc) 





2 


35539E+01 


1 


06202E-31 


1 


35574E- 


-05 


1 


4 


28137E+01 


-1 


05730E-33 


-5 


66429E- 


-07 


2 


6 


81625E+01 


1 


75863E-30 


4 


32438E- 


-07 


3 


9 


96592E+01 


-2 


41634E-32 


1 


02708E- 


-04 


4 


1 


37339E+02 


1 


69583E-25 


-1 


79228E- 


-05 


5 


1 


81223E+02 


3 


64142E-26 


6 


92222E- 


-06 



Table 2: Branching points of the six branches in the Gaussian well example. 




Im(k) lm(k) 



(a) Complex k plane (b) S(fc) X A plane 

Figure 3: Projections of trajectories of the S-matrix poles representing the first six 
bound/resonant states for (I = 3)-waves in a Gauss potential. 
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Figure 4: Schematic view of the full trajectories of the S-matrix poles representing the first 
three bound/resonant states for (I = 3)-waves in a Gauss potential shown in the A; X A space 
along with projections on the various planes. 

4-3. Square potential well 

As a second example we use a s-wave (I = 0) in a square potential well: 



For our purposes we choose a = 1. 

Analytical results for such potentials are well known and were extensively 
studied in [29]. We use them as reference for our numerical studies. 

The grid used for the renormalized Numerov solver was r G [0, 1.1] with 
2048 points and the step size As was the same as in the previous example. 
The first three bound states were used for the continuation and the starting 
points of these three branches are given in table 3. As in the Gauss potential 
case, the continuation was performed from these points in the direction of the 
origin resulting in branches shown in figures 5(a), 5(b) and 6. As follows from 
theoretical considerations in [29], the ground state of this potential does not 
transform into a resonance. All other states n > do branch off into resonances 
at k = —i. The real part of k tends to ±nir as Q(k) — > — oo and A —¥ 0, which 
again, corresponds to theoretical results. 

5. Discussion and Conclusions 

Changing system parameters in the potential of a Schrodinger equation can 
turn resonances into bound states or vice versa. This transition is usually marked 




r > a 



r < a 



(41) 
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n A 



5 2.15040e+00 

1 15 2.02173e+00 

2 32 8.25262e-01 

Table 3: Starting points for the three continuation branches in the square potential well 
example, obtained with the renormalized Numerov method. 5R(fc) = for all points. 




Figure 5: Projections of trajectories of the S-matrix poles representing the first three 
bound/resonant states for s-waves in a square potential well. Note the theoretically well-known 
lack of bifurcation in the ground-state branch. The other bifurcations are located at k = — i 
and not at the origin which is a known result for these types of s-wave problems. 
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Figure 6: Schematic view of the full trajectories of the S-matrix poles for s-waves in a square 
potential well shown in the k X A space along with projections. 

with a double pole of the S-matrix. We have interpreted this double pole as a 
bifurcation point where different states meet. 

In this article we have reviewed the key numerical and mathematical meth- 
ods that allow us to track these poles and detect the bifurcation point in an 
automatic way. These methods, originally developed by the dynamical systems 
community, are based on predictor-corrector methods and are directly applicable 
to the problem at hand. What is needed is a numerical routine that solves the 
radial Schrodinger equation for a given complex k and system parameter A. The 
continuation method then calls this routine multiple times with appropriate ar- 
guments k and A and constructs with this information the complete continuation 
curve. 

Our contribution is the insight to apply the method to a function proportional 
to the Jost function instead of the numerical S-matrix. The latter suffers from 
nearby zeros and poles what can lead to diverging and deteriorating behavior. 
Another insight is to apply the method in the fc-plane rather than the 2?-plane. 
This gives orthogonal tangent directions in the bifurcation point while in the 
_B-plane these are aligned and are much harder to treat numerically. 

Our approach is quite robust. In addition to the examples provided, we 
have applied it to other short-range problems including Morse, Yukawa and 
Lennard- Jones potentials with / ranging from up to 5. In all cases the program 
worked without any modifications. Also for s-wave problems with barriers, where 
the bifurcation does not happen at the origin, the program was able to detect 
the bifurcation point, the tangent directions and follow the solution curves that 
emerge from the bifurcation point. For all these problems the complete curve 
was found within seconds. 
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However, we have found it hard to trace resonances when they move far down 
in the complex plane into the region with large negative imaginary momenta. In 
this region, the fundamental solutions ht and hy are increasing and decreasing 
exponentials with a slight oscillation. In our shooting method that integrates 
outwards starting from r = the exponentially increasing function will dominate 
over the decreasing function. This effect becomes more difficult to deal with 
if R, the end of the domain, is increased. A possible solution to this problem 
might be to use a mismatch function where two shootings, one from the left and 
one from the right, are matched. The shooting from the right would then have 
ht(kR) as boundary condition. It is the question, however, if the mismatch 
function is suitable for pseudo-arclength continuation near the bifurcation point. 

The proposed method is, however, independent of this solver and can be 
built around any solver of the radial Schrodinger equation. In addition to the 
renormalized Numerov solver, we have tested the method with a finite difference 
matrix method with very similar results. 

In the future, we will extend the method to problems with unknown asymp- 
totic solutions that require absorbing boundary conditions such as ECS [23]. 
This will allow us to track resonances in multidimensional scattering problems. 
Note, that it is then not possible to solve for the 5-matrix in the complete 
complex fc-plane. It is for these large scale problems that our method can prove 
to be valuable since other methods, based on interpreting a resonance as an 
eigenvalue, then become intractable because they do not scale to a large number 
of unknowns. Finally, as indicated in section 3 the application of our methods 
to long-range potentials is important as well. This would allow us to tackle 
problems with a broader range of physical applications and to verify our results 
with experimental data. 
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